Visualize geographic map¶
Created by
¶ Li, Chaonan (李超男) / licn@mtc.edu.cn / Ecological Security and Protection Key Laboratory of Sichuan Province, Mianyang Normal University (绵阳师范学院生态安全与保护四川省重点实验室)
¶ Liu Chi (刘驰) / liuchi0426@126.com / College of Resources and Environment, Fujian Agriculture and Forestry University (福建农林大学资源与环境学院)
¶ Liao, Haijun (廖海君) / lihj@mtc.edu.cn / Engineering Research Center of Chuanxibei RHS Construction at Mianyang Normal University of Sichuan Province (绵阳师范学院川西北乡村人居环境建设工程研究中心)
Reviewed by
¶ Li, Xiangzhen (李香真) / lixz@fafu.edu.cn / College of Resources and Environment, Fujian Agriculture and Forestry University (福建农林大学资源与环境学院)
One of our purposes in developing the microgeo R package is to rapidly visualize the biogeographic patterns of microbial traits onto maps. Hence, we have designed three visualization functions, and all of them return a ggplot2 object. The function of plot_bmap() is able to accept a microgeo-compatible SpatialPolygonsDataFrame returned by read_aliyun_map(), trans_map_fmt(), grid_map(), merge_dfs_to_map() or merge_mtx_to_map(), and it only visualize the data included in the SpatialPolygonsDataFrame itself. The plot_nmap() and plot_imap() are designed to visualize the results of nearest neighbour and inverse distance weighting interpolations, respectively.
To add more microbial traits onto maps, we have implemented four functions including add_label(), add_point(), add_spatraster() and add_contour. The add_label() can add any characters or numbers in the native polygons of map or the grids created by grid_map(). For example, we can apply the gradient colors to represent the average shannon index in each grid, and simultaneously add the standard errors into each grid by using the add_label(). The add_point() can add points based on a numeric vector onto a map. By applying such a function, we also can use the size of point to represent the standard errors of shannon index. Meanwhile, we also can use the gradient colors and/or size of point as a proxy of any microbial trait at each sampling site on maps.
In the microgeo R package, the functions responsible for predictions based on machine learning model and some spatial interpolations return a SpatRaster object. To visualize such a type of data onto maps, we have designed a function of add_spatraster(). To convert the data in the SpatRaster to contours that can be further rendered on maps, we have designed a function of add_contours().
Now, let's go through each of these functions and see how to visualize microbial/environmental traits onto maps.
Load required R packages¶
# Load three required packages
suppressMessages(require("magrittr"))
require("ggplot2") %>% suppressMessages()
require("microgeo") %>% suppressMessages()
prev_locale <- Sys.setlocale("LC_CTYPE", "C.UTF-8") # Set locale to `UTF-8` to show Chinese characters in tables
Plot basic map¶
The microgeo R package provides two test datasets for the plot_bmap() function, where the common.map.mean is a SpatialPolygonsDataFrame based on administrative divisions and the gridded.map.mean is a SpatialPolygonsDataFrame processed by grid_map(). Below is a simple example of the plot_bmap() function.
Firstly, Let's check the content of common.map.mean, a microgeo-compatible SpatialPolygonsDataFrame.
# This is a microgeo-compatible SpatialPolygonsDataFrame
data(common.map.mean)
head(common.map.mean@data[,1:12])
| TYPE | FMTS | NAME | X.CENTER | Y.CENTER | observed_mean | shannon_mean | observed_sd | shannon_sd | observed_se | shannon_se | sample.num | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | |
| 1 | DataV.GeoAtlas | microgeo | 西藏自治区 | 88.38828 | 31.56375 | 663.8416 | 5.846949 | 239.0407 | 0.4781792 | 10.44254 | 0.02088936 | 524 |
| 2 | DataV.GeoAtlas | microgeo | 青海省 | 96.04353 | 35.72640 | 647.7036 | 5.837889 | 246.2542 | 0.5396712 | 11.71312 | 0.02566954 | 442 |
| 3 | DataV.GeoAtlas | microgeo | 四川省 | 102.69345 | 30.67454 | 706.8417 | 5.974116 | 198.0067 | 0.3659859 | 11.87565 | 0.02195037 | 278 |
# Change the name of Polygons
common.map.mean$NAME <- c("Tibet", "Qinghai", "Sichuan")
head(common.map.mean@data[,1:12]) %>% suppressMessages()
Loading required package: sp
| TYPE | FMTS | NAME | X.CENTER | Y.CENTER | observed_mean | shannon_mean | observed_sd | shannon_sd | observed_se | shannon_se | sample.num | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | |
| 1 | DataV.GeoAtlas | microgeo | Tibet | 88.38828 | 31.56375 | 663.8416 | 5.846949 | 239.0407 | 0.4781792 | 10.44254 | 0.02088936 | 524 |
| 2 | DataV.GeoAtlas | microgeo | Qinghai | 96.04353 | 35.72640 | 647.7036 | 5.837889 | 246.2542 | 0.5396712 | 11.71312 | 0.02566954 | 442 |
| 3 | DataV.GeoAtlas | microgeo | Sichuan | 102.69345 | 30.67454 | 706.8417 | 5.974116 | 198.0067 | 0.3659859 | 11.87565 | 0.02195037 | 278 |
Then, we use plot_bmap to visualize the common.map.mean.
# Plot a basic map without any fill color
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p1 <- common.map.mean %>% plot_bmap(legend.position = c(0.06, 0.25)) %>% suppressMessages())
# Plot a basic map with the naive polygons filled by the automatically generated colors based on a character variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p2 <- common.map.mean %>% plot_bmap(var = 'NAME', fill = 'auto', legend.position = c(0.06, 0.25)))
# Plot a basic map with the naive polygons filled by the automatically generated colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p3 <- common.map.mean %>% plot_bmap(var = 'shannon_mean', fill = 'auto', legend.position = c(0.06, 0.255)))
# Plot a basic map with the naive polygons filled by manually set colors based on a character variable
# The display order of polygon names are also manually set
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p4 <- common.map.mean %>% plot_bmap(var = 'NAME', fill = RColorBrewer::brewer.pal(8, "Set2")[1:3],
ord = c('Tibet', 'Qinghai', 'Sichuan'), legend.position = c(0.06, 0.25)))
# Plot a basic map with the naive polygons filled by manually set colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p5 <- common.map.mean %>% plot_bmap(var = 'shannon_mean', fill = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100),
legend.position = c(0.06, 0.255)))
Now, let's check the content of gridded.map.mean, a microgeo-compatible SpatialPolygonsDataFrame. The NA means that there are no samples in the grid. The NAME means the grid ID.
# This is a microgeo-compatible SpatialPolygonsDataFrame
data(gridded.map.mean)
head(gridded.map.mean@data[,1:12])
| TYPE | FMTS | NAME | X.CENTER | Y.CENTER | observed_mean | shannon_mean | observed_sd | shannon_sd | observed_se | shannon_se | sample.num | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | |
| 1 | Gridded.Map | microgeo | 1 | 92.72708 | 27.05839 | NA | NA | NA | NA | NA | NA | NA |
| 2 | Gridded.Map | microgeo | 2 | 93.60240 | 27.08612 | NA | NA | NA | NA | NA | NA | NA |
| 3 | Gridded.Map | microgeo | 3 | 102.68131 | 26.73440 | NA | NA | NA | NA | NA | NA | NA |
| 4 | Gridded.Map | microgeo | 4 | 101.88706 | 26.72948 | NA | NA | NA | NA | NA | NA | NA |
| 5 | Gridded.Map | microgeo | 5 | 94.08106 | 28.05720 | NA | NA | NA | NA | NA | NA | NA |
| 6 | Gridded.Map | microgeo | 6 | 100.48112 | 28.28070 | NA | NA | NA | NA | NA | NA | NA |
Then, we use plot_bmap to visualize the gridded.map.mean.
# Plot a gridded map without any fill color
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p6 <- gridded.map.mean %>% plot_bmap(legend.position = c(0.06, 0.25))) %>% suppressMessages()
# Plot a gridded map with the grids filled by the automatically generated colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p7 <- gridded.map.mean %>% plot_bmap(var = 'shannon_mean', fill = 'auto', legend.position = c(0.06, 0.255)))
# Plot a gridded map with the grids filled by manually set colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p8 <- gridded.map.mean %>% plot_bmap(var = 'shannon_mean', fill = rev(colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100)),
legend.position = c(0.06, 0.255)))
Plot a nearest neighbour interpolation map¶
We have designed a function of plot_nmap() to visualize the results of nearest neighbour interpolation. Here, we only give a simple example for the visualization. See details in the tutorial of nearest neighbour interpolation.
Firstly, let's check the content of nen.interp.
# Load test data
data(nen.interp)
# The `nen.interp` is a list
class(nen.interp)
# There are two datasets in the `nen.interp`
names(nen.interp)
- 'shannon'
- 'observed'
Then, we use plot_nmap to visualize the nen.interp.
# Plot a map for the results of nearest neighbour interpolation [alpha diversity ==> observed]
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p9 <- nen.interp$observed %>% plot_nmap(legend.position = c(0.06, 0.25)))
# Plot a map for the results of nearest neighbour interpolation [alpha diversity ==> shannon]; Use manually set colors
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p10 <- nen.interp$shannon %>% plot_nmap(fill = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100), legend.position = c(0.06, 0.255)))
Plot an inverse distance weighting (IDW) interpolation map¶
The microgeo provides two types of inverse distance weighting (IDW) interpolation methods, one of which is hexagonal interpolation. It generates a sf data.frame and requires a separate plotting function for visualization. To address this issue, we have designed a function of plot_imap(). Here, we only give a simple example for visualization. See details in the tutorial of IDW interpolation.
Firstly, let's check the content of idw.interp.hex.
# Load test data
data(idw.interp.hex)
# The `idw.interp.hex` is a list
class(idw.interp.hex)
# There are two datasets in the `idw.interp.hex`
names(idw.interp.hex)
- 'observed'
- 'shannon'
Then, we use plot_imap to visualize the idw.interp.hex.
# Plot a map for the results of inverse distance weighting (IDW) interpolation [alpha diversity ==> observed]
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p11 <- idw.interp.hex$observed %>% plot_imap(legend.position = c(0.06, 0.255)))
# Plot a map for the results of inverse distance weighting (IDW) interpolation [alpha diversity ==> observed]; Use manually set colors
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p12 <- idw.interp.hex$observed %>% plot_imap(fill = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100), legend.position = c(0.06, 0.255)))
Add labels onto a map¶
If you want to add some labels onto maps, the add_label() function may be very useful. With the add_label(), you can add the names of administrative region, grid names, and any character or numeric variables included in SpatialPolygonsDataFrame onto maps. You can even use this function multiple times to add multiple labels onto maps. Here are a few simple examples.
# Add grid NAME on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p13 <- p6 %>% add_label(dat = gridded.map.mean@data, lab.var = 'NAME', lon.var = 'X.CENTER',
lat.var = 'Y.CENTER', size = 5, remove.na = FALSE) # show all grid names
p14 <- p6 %>% add_label(dat = gridded.map.mean@data, lab.var = 'NAME', lon.var = 'X.CENTER',
lat.var = 'Y.CENTER', size = 5, remove.na = TRUE) # only show the names of grid covered by field sampling
cowplot::plot_grid(p13, p14, align = 'hv', ncol = 2, label_size = 28, labels = c("a", "b")) %>% suppressWarnings()
# Add the names of naive polygons and sample number in each naive polygon onto the map
options(repr.plot.width = 15.4, repr.plot.height = 8.02 * 3)
p15 <- p4 %>%
add_label(dat = common.map.mean@data, lab.var = 'NAME', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5) %>%
add_label(dat = common.map.mean@data, lab.var = 'sample.num', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5, vjust = 3) %>% suppressWarnings()
# Add the sample number of each grid onto the map
p16 <- p7 %>% add_label(dat = gridded.map.mean@data, lab.var = 'sample.num', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5, remove.na = F)
# Add the standard errors of shannon index in each grid on the map
# Grids without labels mean that the samples are not enough for the calculation of standard error (at least two samples are required)
p17 <- p7 %>% add_label(dat = gridded.map.mean@data, lab.var = 'shannon_se', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5, remove.na = F)
# Display the figure
cowplot::plot_grid(p15, p16, p17, align = 'hv', ncol = 1, label_size = 28, labels = c("a", "b", "c")) %>% suppressWarnings()
Add points onto a map¶
In some cases, we want to display additional data onto the map processed by grid_map(). For example, we need to display both the mean and standard error of the Shannon index in each grid at the same time. The first option is to represent the mean by the color of the grid, and then use the add_label() function to directly add the value of the standard error to the grid. Another option is to use add_point() function, which adds points to each grid, with the size of the point representing the standard error.
# Add points onto a gridded map by using the standard error of shannon index in each grid
# Grids without point mean that the samples are not enough for the calculation of standard error (at least two samples are required)
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p18 <- p7 %>% add_point(dat = gridded.map.mean@data, lab.var = 'shannon_se',
lon.var = "X.CENTER", lat.var = "Y.CENTER", remove.na = TRUE) )# Remove all NA
You also can use the add_point() to visualize a numeric variable onto maps. For example, we can visualize the soil pH of each sampling site onto maps. The below is a simple example.
# Create a microgeo dataset. See subsequent tutorial for details.
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude")
dataset.dts %<>% tidy_dataset()
ℹ [2024-01-04 17:10:15.457001] INFO ==> all samples fall within the map area! ℹ [2024-01-04 17:10:15.464771] INFO ==> dataset has been created successfully! ℹ [2024-01-04 17:10:15.470771] INFO ==> use `object %>% show_dataset()` to check the summary of dataset.
# Visualize the soil pH of each sampling site onto maps
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02 * 3)
dat <- cbind(dataset.dts$met, dataset.dts$env)
p19 <- p1 %>% add_point(dat = dat, lab.var = 'pH')
p20 <- p1 %>% add_point(dat = dat, lab.var = 'pH', fill = 'lightblue')
p21 <- p1 %>% add_point(dat = dat, lab.var = 'pH', fill = rev(colorRampPalette(RColorBrewer::brewer.pal(11, 'RdYlGn'))(100)))
p22 <- p7 %>% add_point(dat = dat, lab.var = 'pH')
p23 <- p9 %>% add_point(dat = dat, lab.var = 'pH')
p24 <- p11 %>% add_point(dat = dat, lab.var = 'pH')
cowplot::plot_grid(p19, p20, p21, p22, p23, p24, align = 'hv', ncol = 2, label_size = 28,
labels = c("a", "b", "c", "d", "e", "f")) %>% suppressWarnings()
Add SpatRasters onto a map¶
In some cases, we need to visualize SpatRaster data. For example, spatial datasets (such as aridity index, annual average temperature, and soil pH) downloaded from public databases, and the results returned by some spatial interpolation functions and machine learning prediction functions in the microgeo R package. Regarding to the visualization of SpatRaster data, we have designed a function of add_spatraster(). Here are a few simple examples.
# Create a microgeo dataset. See subsequent tutorial for details.
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude")
ℹ [2024-01-04 17:10:23.413691] INFO ==> all samples fall within the map area! ℹ [2024-01-04 17:10:23.421172] INFO ==> dataset has been created successfully! ℹ [2024-01-04 17:10:23.427492] INFO ==> use `object %>% show_dataset()` to check the summary of dataset.
# Download phh2o from soilgrid database.
dataset.dts %<>% get_soilgrid(depth = 5, measures = c('phh2o'), out.dir = "../dev/dat/rundata")
✔ [2024-01-04 17:10:24.499287] SAVE ==> results have been saved to: object$spa$rast$his$phh2o
# Add the SpatRaster data of soil pH onto a map
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p25 <- p1 %>% add_spatraster(spat.raster = dataset.dts$spa$rast$his$phh2o))
Notably, the add_spatraster() function would fill the whole study area. However, in some cases, we only focus one or two ecosystems, for example, grassland, forest or barren. To address this problem, we also have implemented a function of mask_spatraster_by_cla(), which can mask off those areas that we are not interested in. Here are a few simple examples.
# Download classification MODIS metrics of research region.
# Please replace the <username> and <passwd> with real value
dataset.dts %<>% get_modis_cla_metrics(username = "username", password = "passwd", out.dir = "../dev/dat/rundata")
ℹ [2024-01-04 17:10:27.520267] INFO ==> preparing MODIS product list for searching... ℹ [2024-01-04 17:10:27.535475] INFO ==> searching avaliable MODIS products... ℹ [2024-01-04 17:10:27.542682] INFO ==> current product (1/1): MCD12Q1 (LC_Type1--> 2022-01-01 to 2022-12-31) ℹ [2024-01-04 17:10:29.238166] INFO ==> find 8 files with 0.09 GB in total... ℹ [2024-01-04 17:10:29.247663] INFO ==> downloading all avaliable MODIS products[skip if the file exists]... ℹ [2024-01-04 17:10:29.257661] INFO ==> preparing the PTVs (Product, Time, Version) for merging remote-sensing images... ℹ [2024-01-04 17:10:29.28531] INFO ==> converting hdf files to tif files... ℹ [2024-01-04 17:10:29.293728] INFO ==> current product (1/1): MCD12Q1 (convert 8 hdf files into 1 tif files using 1 threads) ℹ [2024-01-04 17:10:29.506624] INFO ==> collecting all merged image files... ℹ [2024-01-04 17:10:29.515717] INFO ==> current measure (1/1): LC_Type1_061 ℹ [2024-01-04 17:10:29.588037] INFO ==> reprojecting the CRS of SpatRaster to epsg:+proj=longlat +datum=WGS84 +no_defs, it takes a while... ✔ [2024-01-04 17:11:01.754878] SAVE ==> results have been saved to: object$spa$rast$cla$LC_Type1
# Mask the results by using grassland(10); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m1 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = 10)
(p26 <- p1 %>% add_spatraster(spat.raster = m1))
# Mask the results by using barren(16); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m2 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = 16)
(p27 <- p1 %>% add_spatraster(spat.raster = m2))
# Mask the results by using forest (1,2,3,4,5); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m3 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = c(1,2,3,4,5))
(p28 <- p1 %>% add_spatraster(spat.raster = m3))
# Mask the results by using grassland (10) and barren (16); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m4 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = c(10, 16))
(p29 <- p1 %>% add_spatraster(spat.raster = m4))
Add contours onto a map¶
You also can convert the SpatRaster returned by spatial interpolation functions or machine learning prediction functions to contours, and further add these contours onto maps by using add_contour(). Here are a few simple examples.
# Create a microgeo dataset. See subsequent tutorial for details.
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude")
dataset.dts %<>% rarefy_count_table()
dataset.dts %<>% tidy_dataset()
dataset.dts %<>% calc_alpha_div(measures = c("observed", "shannon"))
ℹ [2024-01-04 17:11:48.354824] INFO ==> all samples fall within the map area! ℹ [2024-01-04 17:11:48.364121] INFO ==> dataset has been created successfully! ℹ [2024-01-04 17:11:48.372576] INFO ==> use `object %>% show_dataset()` to check the summary of dataset. ℹ [2024-01-04 17:11:50.699992] INFO ==> the ASV/gene abundance table has been rarefied with a sub-sample depth of 5310 ✔ [2024-01-04 17:11:53.637593] SAVE ==> results have been saved to: object$div$alpha
# Perform kriging interpolation for shannon diversity indices.
kri.rst.shannon <- interp_kri(map = dataset.dts$map, met = dataset.dts$met, dat = dataset.dts$div$alpha,
var = 'shannon', model = 'Mat', trim.dup = TRUE)
! [2024-01-04 17:11:53.963551] WARN ==> only use 420 out of 1244 sampling sites for interpolation!
[using ordinary kriging]
# Add contours on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p30 <- p1 %>% add_contour(spat.raster = kri.rst.shannon, nlevels = 100)
p31 <- p1 %>% add_spatraster(spat.raster = kri.rst.shannon) %>% add_contour(spat.raster = kri.rst.shannon, nlevels = 10, color = 'white', show.labels = TRUE)
cowplot::plot_grid(p30, p31, align = 'hv', ncol = 2, label_size = 28, labels = c("a", "b")) %>% suppressWarnings()
Add north arrow onto a map¶
We have designed a function of add_north_arrow() to add a north arrow onto the map. You can use ?add_north_arrow() to view the built-in document. Here are a few simple examples.
# Add a north arrow on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p32 <- p29 %>% add_north_arrow() %>% suppressWarnings()
p33 <- p30 %>% add_north_arrow() %>% suppressWarnings()
cowplot::plot_grid(p32, p33, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings()
Add spatial-aware scale bar onto a map¶
We have designed a function of add_scale_bar() to add a spatial-aware scale bar onto the map. You can use ?add_scale_bar() to view the built-in document. Here are a few simple examples.
# Add a spatial-aware scale bar on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p34 <- p32 %>% add_scale_bar()
p35 <- p33 %>% add_scale_bar()
cowplot::plot_grid(p34, p35, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings() %>% suppressMessages()
Add CRS onto a map before final visualization¶
We have designed a function of add_crs() to add a coordinate reference system to the map before final visualization. You can use ?add_crs() to view the built-in document. Here are a few simple examples.
# Add a coordinate reference system to the map before final visualization
options(repr.plot.width = 13.4 * 2, repr.plot.height = 8.02)
p36 <- p34 %>% add_crs()
p37 <- p35 %>% add_crs()
cowplot::plot_grid(p36, p37, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings() %>% suppressMessages()
Add sampling sites onto a map¶
In the large-scale microbial biogeographical studies, another application of microgeo R package is to draw sample point distribution maps. In the current version, we provide a flexible function of add_sampl_site(), which can render sampling points on the map. Here are a few simple examples.
# Create a microgeo dataset. See subsequent tutorial for more details.
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude")
dataset.dts %<>% get_ai(out.dir = "../dev/dat/rundata")
dataset.dts %<>% get_his_bioc(res = 10, out.dir = "../dev/dat/rundata")
dataset.dts %<>% get_soilgrid(depth = 5, measures = c('phh2o'), out.dir = "../dev/dat/rundata")
dataset.dts %<>% extract_data_from_spatraster(type = 'his')
dataset.dts %<>% tidy_dataset()
ℹ [2024-01-04 17:12:08.654725] INFO ==> all samples fall within the map area! ℹ [2024-01-04 17:12:08.664679] INFO ==> dataset has been created successfully! ℹ [2024-01-04 17:12:08.674071] INFO ==> use `object %>% show_dataset()` to check the summary of dataset. ✔ [2024-01-04 17:12:12.710911] SAVE ==> results have been saved to: object$spa$rast$his$AI ✔ [2024-01-04 17:12:20.137347] SAVE ==> results have been saved to: object$spa$rast$his(19 variables) ✔ [2024-01-04 17:12:26.392774] SAVE ==> results have been saved to: object$spa$rast$his$phh2o ✔ [2024-01-04 17:12:26.411119] SAVE ==> results have been saved to: object$spa$tabs
# Add two grouping variable in `dataset.dts$met`
dataset.dts$met$GA <- cut(dataset.dts$spa$tabs$AI, breaks = c(-Inf, 0.2, 0.5, Inf), labels = c("Arid", "Semi-arid", "Humid"))
dataset.dts$met$GM <- cut(dataset.dts$spa$tabs$Bio1, breaks = c(-Inf, -1, 1, Inf), labels = c("Cold", "Moderate", "Warm"))
# Add sampling sites without colors
options(repr.plot.width = 13.4 * 2, repr.plot.height = 8.02 * 3)
s1 <- p4 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s2 <- p7 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s3 <- p11 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s4 <- p9 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s5 <- p25 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s6 <- p29 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>%
suppressWarnings() %>% add_scale_bar() %>% add_north_arrow() %>% add_crs()
cowplot::plot_grid(s1, s2, s3, s4, s5, s6, align = 'hv', ncol = 2, labels = c("a", "b", "c", "d", "e", "f"), label_size = 28) %>% suppressWarnings()
# Add sampling sites with color and shapes
options(repr.plot.width = 13.4 * 2, repr.plot.height = 8.02)
s7 <- p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", color.val = "gray30",point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s8 <- p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GA", color.val = "gray30", point.alpha = 0.5) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
cowplot::plot_grid(s7, s8, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings()
# Facet by one grouping variable
options(repr.plot.width = 13.43 * 2, repr.plot.height = 7.9 * 3)
s9 <- p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GA", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.col.nums = 1) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
s10 <- p29 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GA", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.col.nums = 1) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
cowplot::plot_grid(s9, s10, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings()
# Facet by two grouping variables [on basic map]
options(repr.plot.width = 13.43 * 3, repr.plot.height = 7.23 * 3)
p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GM", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.by.shape = T) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()
# Facet by two grouping variables [on SpatRaster]
options(repr.plot.width = 13.43 * 3, repr.plot.height = 7.9 * 3)
p29 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GM", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.by.shape = T) %>%
add_scale_bar() %>% add_north_arrow() %>% add_crs()